function [TT_step2, TT_deriv_1_step2, TT_deriv_2_step2,TT_deriv_3_step2] = construct_chebyshev_berkson_Ponly_twodim_inc(x_interest_1,x_interest_2,x_interest_3,K1,K2,K3,a1,a2,a3,b1,b2,b3,trunc_low1,trunc_low2,trunc_low3,trunc_high1,trunc_high2,trunc_high3,mu1,mu2,mu3,sigma1,sigma2,sigma3,fdist1,fdist2,fdist3,alpha1,alpha2,alpha3)

[T_1, T_deriv_1]=construct_chebyshev_berkson(x_interest_1,K1,a1,b1,trunc_low1,trunc_high1,mu1,sigma1,fdist1,alpha1);
[T_2, T_deriv_2]=construct_chebyshev(x_interest_2,K2,a2,b2);
[T_3, T_deriv_3]=construct_chebyshev(x_interest_3,K3,a3,b3);

TT=[];
TT_deriv_1=[];
TT_deriv_2=[];

for ii=1:(K2+1)

    T_add=T_1.*(T_2(:,ii)*ones(1,1+K1));
    TT=[TT T_add];    

    T_deriv_add1=T_deriv_1 .* (T_2(:,ii)*ones(1,1+K1));
    TT_deriv_1=[TT_deriv_1 T_deriv_add1];
    
    T_deriv_add2=T_1 .* (T_deriv_2(:,ii)*ones(1,1+K1));
    TT_deriv_2=[TT_deriv_2 T_deriv_add2];
    
end

TT_step2=[];
TT_deriv_1_step2=[];
TT_deriv_2_step2=[];
TT_deriv_3_step2=[];

for ii=1:(K3+1)

    T_add=TT.*(T_3(:,ii)*ones(1,(1+K1)*(1+K2)));
    TT_step2=[TT_step2 T_add];    

    T_deriv_add1=TT_deriv_1 .* (T_3(:,ii)*ones(1,(1+K1)*(1+K2)));
    TT_deriv_1_step2=[TT_deriv_1_step2 T_deriv_add1];
    
    T_deriv_add2=TT_deriv_2 .* (T_3(:,ii)*ones(1,(1+K1)*(1+K2)));
    TT_deriv_2_step2=[TT_deriv_2_step2 T_deriv_add2];
    
    T_deriv_add3=TT .* (T_deriv_3(:,ii)*ones(1,(1+K1)*(1+K2)));
    TT_deriv_3_step2=[TT_deriv_3_step2 T_deriv_add3];
    
end

    %
